{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scalable Multitask GP Regression (w/ KISS-GP)\n",
    "\n",
    "This notebook demonstrates how to perform a scalable multitask regression.\n",
    "\n",
    "It does everything that the [standard multitask GP example](./Multitask_GP_Regression.ipynb) does, but using the SKI scalable GP approximation. This can be used on much larger datasets (up to 100,000+ data points).\n",
    "\n",
    "For more information on SKI, check out the [scalable regression example](../04_Scalable_GP_Regression_1D/KISSGP_Regression_1D.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up training data\n",
    "\n",
    "In the next cell, we set up the training data for this example. We'll be using 1000 regularly spaced points on [0,1] which we evaluate the function on and add Gaussian noise to get the training labels.\n",
    "\n",
    "We'll have two functions - a sine function (y1) and a cosine function (y2).\n",
    "\n",
    "For MTGPs, our `train_targets` will actually have two dimensions: with the second dimension corresponding to the different tasks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_x = torch.linspace(0, 1, 1000)\n",
    "\n",
    "train_y = torch.stack([\n",
    "    torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * 0.2,\n",
    "    torch.cos(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * 0.2,\n",
    "], -1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the model\n",
    "\n",
    "The model should be somewhat similar to the `ExactGP` model in the [simple regression example](../01_Simple_GP_Regression/Simple_GP_Regression.ipynb).\n",
    "\n",
    "The differences:\n",
    "\n",
    "1. We're going to wrap ConstantMean with a `MultitaskMean`. This makes sure we have a mean function for each task.\n",
    "2. Rather than just using a RBFKernel, we're using that in conjunction with a `MultitaskKernel`. This gives us the covariance function described in the introduction.\n",
    "3. We're using a `MultitaskMultivariateNormal` and `MultitaskGaussianLikelihood`. This allows us to deal with the predictions/outputs in a nice way. For example, when we call MultitaskMultivariateNormal.mean, we get a `n x num_tasks` matrix back.\n",
    "\n",
    "In addition, we're going to wrap the RBFKernel in a `GridInterpolationKernel`.\n",
    "This approximates the kernel using SKI, which makes GP regression very scalable.\n",
    "\n",
    "You may also notice that we don't use a ScaleKernel, since the IndexKernel will do some scaling for us. (This way we're not overparameterizing the kernel.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MultitaskGPModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(MultitaskGPModel, self).__init__(train_x, train_y, likelihood)\n",
    "        \n",
    "        # SKI requires a grid size hyperparameter. This util can help with that\n",
    "        grid_size = gpytorch.utils.grid.choose_grid_size(train_x)\n",
    "        \n",
    "        self.mean_module = gpytorch.means.MultitaskMean(\n",
    "            gpytorch.means.ConstantMean(), num_tasks=2\n",
    "        )\n",
    "        self.covar_module = gpytorch.kernels.MultitaskKernel(\n",
    "            gpytorch.kernels.GridInterpolationKernel(\n",
    "                gpytorch.kernels.RBFKernel(), grid_size=grid_size, num_dims=1,\n",
    "            ), num_tasks=2, rank=1\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultitaskMultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "    \n",
    "likelihood = gpytorch.likelihoods.MultitaskGaussianLikelihood(num_tasks=2)\n",
    "model = MultitaskGPModel(train_x, train_y, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train the model hyperparameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/50 - Loss: 404.388\n",
      "Iter 2/50 - Loss: 352.338\n",
      "Iter 3/50 - Loss: 300.768\n",
      "Iter 4/50 - Loss: 249.549\n",
      "Iter 5/50 - Loss: 198.132\n",
      "Iter 6/50 - Loss: 145.569\n",
      "Iter 7/50 - Loss: 91.422\n",
      "Iter 8/50 - Loss: 38.138\n",
      "Iter 9/50 - Loss: -11.845\n",
      "Iter 10/50 - Loss: -61.598\n",
      "Iter 11/50 - Loss: -110.397\n",
      "Iter 12/50 - Loss: -157.206\n",
      "Iter 13/50 - Loss: -203.842\n",
      "Iter 14/50 - Loss: -249.601\n",
      "Iter 15/50 - Loss: -292.015\n",
      "Iter 16/50 - Loss: -340.180\n",
      "Iter 17/50 - Loss: -384.880\n",
      "Iter 18/50 - Loss: -427.049\n",
      "Iter 19/50 - Loss: -471.487\n",
      "Iter 20/50 - Loss: -516.107\n",
      "Iter 21/50 - Loss: -553.546\n",
      "Iter 22/50 - Loss: -600.415\n",
      "Iter 23/50 - Loss: -634.620\n",
      "Iter 24/50 - Loss: -676.436\n",
      "Iter 25/50 - Loss: -715.459\n",
      "Iter 26/50 - Loss: -754.357\n",
      "Iter 27/50 - Loss: -792.451\n",
      "Iter 28/50 - Loss: -826.312\n",
      "Iter 29/50 - Loss: -854.712\n",
      "Iter 30/50 - Loss: -887.396\n",
      "Iter 31/50 - Loss: -918.314\n",
      "Iter 32/50 - Loss: -945.657\n",
      "Iter 33/50 - Loss: -971.053\n",
      "Iter 34/50 - Loss: -992.324\n",
      "Iter 35/50 - Loss: -1014.124\n",
      "Iter 36/50 - Loss: -1035.698\n",
      "Iter 37/50 - Loss: -1052.513\n",
      "Iter 38/50 - Loss: -1065.371\n",
      "Iter 39/50 - Loss: -1074.025\n",
      "Iter 40/50 - Loss: -1083.747\n",
      "Iter 41/50 - Loss: -1090.505\n",
      "Iter 42/50 - Loss: -1092.604\n",
      "Iter 43/50 - Loss: -1097.734\n",
      "Iter 44/50 - Loss: -1096.438\n",
      "Iter 45/50 - Loss: -1098.831\n",
      "Iter 46/50 - Loss: -1098.317\n",
      "Iter 47/50 - Loss: -1094.847\n",
      "Iter 48/50 - Loss: -1092.078\n",
      "Iter 49/50 - Loss: -1088.565\n",
      "Iter 50/50 - Loss: -1082.591\n"
     ]
    }
   ],
   "source": [
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam([\n",
    "    {'params': model.parameters()},  # Includes GaussianLikelihood parameters\n",
    "], lr=0.1)\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "n_iter = 50\n",
    "for i in range(n_iter):\n",
    "    optimizer.zero_grad()\n",
    "    output = model(train_x)\n",
    "    loss = -mll(output, train_y)\n",
    "    loss.backward()\n",
    "    print('Iter %d/%d - Loss: %.3f' % (i + 1, n_iter, loss.item()))\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make predictions with the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAADOCAYAAAB//QjSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzsnXl4VNXdgN87k2UCWYYQsiJkY1WBQBCwohjiBgKCIBawtGXr4lYpVgXaWhS0rW1BbRXDZ61gjQQEVFyJAkEWg1BUAgIJgpIwgZCwZZ/7/XFmbmYmM8kkM0km5LzPkyczd86998yde+7vnN+qqKqKRCKRSCSS1kPX1h2QSCQSiaSjIYWvRCKRSCStjBS+EolEIpG0MlL4SiQSiUTSykjhK5FIJBJJKyOFr0QikUgkrUyHF76KohgVRZmrKMpgRVHSFUWZa/NZuqIoa1u5P88qivKoi8/SFUU559DHRxVFeVlRFKOTth+3QP8mW/4PVhRlr6W/iTafJ1qvmTvXz7afTvb1Sv9tr6nl957sjeNKfA85npvcPzme24gOL3yBtcBbqqp+qarqJ0CJ9Ye1vG9tMl19YOnPSofNX6qqOk9V1VInbUvxIoqipANfWo7/JZAPZKqqmm9z3nxVVafY9KFBbPvpZF9v9V+7ptbrZPuAkVxRyPHsJnI8ty0dWvgqijIY6n5Ay+ss4HGbZomWWdtky81qnSWm2/wZLTPWdMusO11RlI8t+zxq+b/X0i5dUZTXbNtbjvmo5fjpjXT7ZWCezXujZX/t3E6+p+1s9FFFUZ61vHbst933cnLuW2wHpqtr6myGa/n+kx3P2di+lnaPWn8ryzZrX+c2ss3pNbX8xrbXUHIFIMezHM/tiQ4tfIFUxGyvHkqd2qdEVdVPLD/ws5ZtU0GbzeUjBvcnlvdDLP8TVVXNUlX1z5Z9SywPhRKgyLa95QazztQbnF1aB4siVDpGIF8Rs75EVVVXAr9zso/trDPL5iO7fjv5Xo4YnWxzPNeXOMxwLQNmis2D0PacDe2baGmXZe2bIlYxudaZumWQOtvW2DVtdzNlSaPI8SzHc7uhowvffJz/aEab2bPtzVNiGRjLgFsURTmGuIEHA+GW2dzLlrZfOhxzreWmDXfS/hZcPDRcYJ0tp1vUa/mqqq5UHOxEbuDYD8fv1SRczK4TLX21Hs/ZtXKFs2tyC3W/Sb7lvattDV3TkkbOLWl/yPEsx3O7oUMLX8ssKtz2JleE8f7PNs1sb9pSy0w1XVXV3yFmeunAx5bjWe0mzngLcdOWOGn/BeImBjGYG+v3SkBzMrCoZ5w6dTjB9uHk2A/H79VUBjvZlgvMoW6V4c61aogvqfsOiYhr52xbk66ppP0jx7Mcz+0Jv7bugA8wGpirKIp1ZhtuuWGtWNVARupmdkMVRbF+nqWqar7FHgGA5f9gRVEGW25KVFUtVRSlxPL+S9v2qqr+2cYOMhgxW13p6HThgK36JRExU0y09HcyYiAMVhQl0fKA+cIykzUC6YqiGG3Oaz2m3fdyck6tPzZ9nWq5PrcApdbtVtUZQhUYbvnsWVVVf+dwrUps2htdvE63vDY67D9YVdU/W47japura+pV5xWJzyDHsxzP7QJFlVWNJG5iGez5jTlp+DpXyveQSDzhShkH7fV7dGi1s6RpWNR6zlRR7QarSrK9DVSJxNvI8dy2yJWvRCKRSCStjFz5SiQSiUTSynjscJWammr1orslNze3XkyaRCJpH8ixLJG0Hh6tfFNTUwcjBuonwODU1NR2F+gskUjkWJZIWhuv2HxTU1ONwLO5ubkuU3w99thj0rgskbjJM888ozTeyvvIsSyReBdXY9lbcb6puBFn9eSTTzb4uclkIjIy0ktd8g6yT+4h++Qe7vTpD3/4Qyv1xileGcvQfq9/a+Jr/QHZJ3fxdCx7xeHKoqoypqamtsvSThKJRCDHskTSOnhq8302NTXVWnWilCsk7ZdE0tGQY1kiaV08VTu/DCRavCSNubm5jrUpJVcANTU1FBUVUVlZiSsfAbPZzPnz51u5Zw3j631SFIXAwECio6Px82vzTK9yLEskrYhHIz43NzefumTabVGoWtIKFBUV0blzZ+Li4rDJFWtHdXU1/v7+rdyzhvH1PqmqSmlpKUVFRXTv3r1N+yXHskTSusgkG5JGqaysxGg0uhS8kuahKApGo5HKysq27opEImllpPCVNIqqqlLwthCKorhU5UskkisXKXwlXqWwsJD09HSKiorauisSiUTis0jhK/Eqy5Yt4/PPP2fp0qXNPsa+fftYv349W7ZsISMjg/x8YYpcv349Cxcu9FZXXVJaWsqYMWOc9qtfv35s2bKFLVu28Nxzz1Fa2u7KiEokEh9ACl+JV4iIiMBgMLBy5UrMZjMrV67EYDBgNBqbdJzS0lIyMjKYNGkSo0ePZvbs2SxatAiAtLS0luh6PYxGIwkJCfW2p6SkkJCQwOjRoxk9ejTz589n2rRpTo9RWlrKc88919JdlUgk7ZQ2j2+QXBl8/fXXLFq0iE2bNlFeXk5QUBATJkzgmWeeadJxsrKySElJsdvWpUsX9u3bR0JCAvv27WPfvn1kZ2cza9Ys9u7dS3h4ONnZ2UycOJHs7GzCw8MZNGgQe/fu5e233yYhIYE+ffrwwQcf8MYbb/DrX/+a+fPnA9i1Dw8PJysri4SEBAoKCtzqr9FopLS0lJKSErKzsykrK2PixIkUFBSwd+9e9u3bR1hYmPbZuHHj6N27d5OuiUQiufKQwlfiFaKjowkNDaWyshKDwUBlZSWhoaFER0d79TwpKSmacF61ahUAJSUlTJ48mUWLFjFr1izCw8NZtWoVDz/8MJ999hlPP/00AB988AEAkyZNIjExkWnTptm1Ly0tZf78+SQmJrJ+/Xq3+1RSUkJiYiLh4eGawJ81axbZ2dlaX62fbdiwgUcffdSbl0QikbRDpNpZ4jVMJhNz5sxh27ZtzJkzh9OnTzf5GJMnTyY7O9tuW0FBQb3VsJWJEycya9Ys/vKXv1BVVcWQIUNISUnRBK6t2jstLY3nnnuOIUOGaNsc24eHi8RO586dc6u/paWlJCYmsmXLFm0yYEt+fr7LzyQSScdFrnwlXiMzM1N7vXz58mYdw2g0smDBAjIyMjT17wsvvGDXxqp2nj9/Ps899xxpaWlMnjyZQYMGaWrj8PBwiouLKSgooLS0FKPRyOTJk1m4cKEmkJ966im79vPnzycrK4shQ4ZQUFDAvn377IR+fn4+BQUFbNmyBYD9+/drfSsoKCAhIYGysjKtnXW77WfOjiuRSDoeUvhKfA5b1bItRqNRW6FaP7fabq3vExPrytBWV1dz66232u3/4osvau8TExPt2tvuv3PnznrnT0xMJC8vT3s/evRo7fXs2bO115MmTbLrky3jxo3zuaxbEomk9ZFqZ4lEIpFIWhkpfCUSiUQiaWWk8JVIJBKJpJWRwlcikUgkklZGCl+JRCKRSFoZKXwlEolEImllpPCV+BT5+fmMGDGCffv2NbhNIpFI2jMyzlfiNgZDYAOfNvRZfSoqnBeQT0xM1JJsWGNyy8rKSEhIkIkpJBLJFYMUvhKfIywszOVn+fn5dsUQgHoFDZ577jmefvppMjMzm1zYQSKRSFoDj4RvamqqEUi3vB2am5v7O8+7JPFVXK1WQWST8mbmpkmTJrF+/fp6K17H4glPP/20XUGD+fPns2rVKlJSUli5cqXX+nOlI8eyRNK6eGrzvQcIz83NzQJITU2d63mXJBKRujErK8vpZ7bFEJwVLejSpUtrdPFKQ45liaQV8Wjlm5uba7u0SARebqi9yWRq8HilpaWedKdFkH0Cs9lMdXV1g21qa2u9cq79+/fzyiuvMGDAAAYNGkTnzp354osvyM/P54svvuCPf/wjmZmZJCQk0KVLF44ePUqPHj0oKSnh6NGjvP766+Tn52tFFb744gtNPe0LOLtOZrO50bHR0nh7LIMcO+7ga/0B2Sd38bRPXrH5pqamJgIlubm5+Q21i4yMbPRY7rRpbTp6n86fP++WStkbauehQ4fy5ptvAtjVvd21a5f22rYY/dChQ7XXU6ZMAeC+++4D4J133vHJIgaOfdLpdD5zj3lzLDelXWvia33ytf6A7JO7eNInb4UaTc7NzZ3npWNJJJK2Q45liaQV8Fj4pqamTs7Nzf2z5XV6Y+0lEolvIseyRNJ6eOrtnA48m5qa+rhlk/SQlEjaIXIsSySti6cOV58ASV7qi0QiaSPkWJZIWheZXlIikUgkklZGCl+JRCKRSFoZmV5S0iSe//SY0+21tWb0evfmcg/cLLWbEomkYyNXvhKfY9++fWRkZLBlyxbWr1/PwoUL3d43IyODffv21duntLSUMWPGeLurEolE0izkylfiU5SWlvKXv/yFN954Q9u2d+9et/ctKChg9uzZJCQk2H1mNBrrbZNIJJK2Qq58HSgsLCQ9PZ2ioqK27kqHJCsri7S0NLttCxYsID8/n/Xr12ur4S1btjBmzBi7Ve7evXspKChgy5YtTJs2DRAC2bqKLigoAERlpIyMDNavX09+fr7TYwE899xz2rkc95FIJBJPkMLXgWXLlvH555+zdOnSJu9rFdwHDhyQAtyLGI1GFi1aRFpamlZwYfTo0RiNRlJSUrQcq0OGDCEhIYHRo0drq9yFCxdq+1m3LVq0iISEBBISEli1apXTY2VkZJCWlkZaWhplZWX19pFIJBJPkMLXgtFoxGAwsHLlSsxmMytXrsRgMNCrV696bV0JWavgnjlzZj0BLlfU7jF58mSys7Pttm3ZssVpW3erF4WHhwNw7tw5bZttZSRnx8rOziYhIQGj0cikSZOc7iNpHHnfSyTOkTZfC3l5eTz22GNs2rSJ8vJyFEVh/PjxLF68uF5bWyF7+PBhkpKS7KrV5OXlAbBy5UpNiP/kJz/RBPKKFSu0toWFhdx333387W9/45FHHmH16tVER0e3/BduJq48lb1Vz9doNLJgwQIyMjJISEigrKyMtLQ0EhISyMrKIiEhgcmTJ7Nv3z4KCgo0O+++ffsoKSmx21ZQUMD8+fPJyspiyJAhWrunnnpKO5ZVMDsea8GCBXZtHPexrTEscY2tJsn2vrfFOgZ8/d6XdCxa+r6UwteGrVu3UlFRgV6vp7a2lm+//VarWlFYWEhiYiKqqmrtrULWFQaDgaqqKioqKrTC7rYCubS0tJ4gb+gh1VFISUmpJ9yMRiOJiYl22zZv3mz3H0QdYOs264TAut/OnTu1du4cy7EPjvtIXGM0GqmoqNDeO973trgjoCWS1qal78sOrXa2VYktW7aMoqIiVFXVVrF5eXn06NGDsLAwRowYgaqqhIaGun38qqoqpk2bxtSpUwkKCgJAURQmTJiAqqp2au68vDw7dbfRaGyR7yyRtAZ5eXlMmDABg8EAiInovffey6FDh7Q2rkw91ntfqqwlbUFj96W36NDCd9myZeTk5BAfH6+tTG0JCgpCp9NRWVmpPQDOnz/f6HFDQ0MZMWIEc+bMobi4mK1bt1JeXg6Aqqp8++23HDp0iKlTp2oPJyvOHlJtjaIodit+ifdQVRVFUdq6G14nJiYGvV6vrX4rKirw8/NjxowZ2ljKy8uzm5gGBQXZ3fueOD9KJM2lsfvSW3RItbOjSswWq8o5ICBAE5iu0Pv5UVtTo73X6XSoqsr4SZP55wsvEOCn48EHH+TDDz+02y8vL4+EhAT0ej2qqmrn1Ov1VFZW8tlnn3n8Hb1JYGAgpaWlGI3GK1JQtBWqqlJaWkpgYGBbd8XrOBtjq1evBrBT42VnZ1NeXk5gYCCVlZWEhobSt29ft1XWEom3iYmJITQ0lMrKSgwGg3Zfetvu2yGFr6NzlRWrEOzXrx89knrz4bsbnR9AUUBVMZvNdpuv+VE6nUPD2f/tCcK7GKmprnLZB51OR3p6OvHx8Rw+fBiTyURkZCQmk4mDBw/6lP0rOjqaoqIiSkpKXK6AzWYzOp1vKVJ8vU+KohAYGHjFORk1NLmFOmFqy5133klERARFRUXk5eXx0EMP8e6772I2mwkKCmLChAk888wzLd11SQfF0fHVz8+Prl27EhUVRUpKCqdPnwbgf//7H7fccguffPKJx+O2Qwpf68zGcWVra+t15UwVFBxK+UWhelYdhO+B7R+xeM1nhIZ3482/Ps7eLZtc9mHixImsWbOGwsJCpk6dSl5eHgcPHtQ+96XZvp+fH927d2+wjXXy4EvIPrUNeXl5PPjgQ7zzjuv735F169YBwuwSExOj+UBYtUEtsfKQSKxYTRzTp0/nyJEjGAwGKioqKC4u5vrrr+dPf/oTN954IwcOHKCiooKZM2fW02g2Fd9aFrQShYWFrF+/vi7doI0qVe8fQNKAYS73tQpeVzx1XxoL7ri6QcEL4mFzoaKGZcuWsWfPHlRVJSkpqcXtDBJJSxMTE8M3zbhvR40apan6jhw5AogJsdlsJiMjw+k+0ilL4g7/+9//iIyM5MCBA9q2wsJCgoKC7JyrrPedo9kjISGBPXv2aNutzriOPjtNoUMIX8cBmpyczNmzZ7V0g9ioUmurqzh2YHezz+W4Gm6IbsZgO/XbsWPHtNW4nO1L2iNWT9H8I982ed/PPvvMpVnDz8+5km7x4sXk5OSwaNGiJp9P0nH46U9/yvnz55k5c6a2bfHixR47ktqGJjaVK174FhYWMmLECHbs2EFiYiIGg8EuIUaboegI6xqFotT/CcaOHcucOXM0OwPIGb6kfWD1FA0MbP6KwBn33HOP3XurkLc6ca1evVqG6EnqYTAYMBgMmhkxLy9P22a9dzxhw4YNzd73iha+1ko21vhdRwepNkU1U3b2NKpav0/vvfceL7/8Mh9++KEmdBcvXizDLiQ+T52nqGuHq+ZgFa6hoaGkp6e7HMsyJE5iy+7du+nRo0eLHd+TGGCPhW9qamp6amrqx54ex9s05nHZHti2bRvJycnk5OSwevVqLeC7R48ecoYv8TreGMtGo5FXXnnFW12qR1xcHJ9//jl33nlnPXtbcnIyhw8fbrFzS9ofAwcO5OTJky16jgkTJjTLN8dj4Zubm/uJp8doCQ4ePMjNt3m7eHoXoA8wHLgDuAcYD4y2bEsG9F4723XXXedURX7XXXdJRyyJ1/HGWG4s5aqnFBQUYDabWb9+vTa5DggIAKCmpkb6SEgAezNdS2tDNm7cSN++fZu8X6uGGplMpgY/91ZIzcnSSr44eYGvvvHkQRABpAPXAdcA1wLuDOwq4BhwGNgDbAH2Ap7bma1JPAICAtDpdI1ez9akrcOhnCH71HK4uve0CIJWpKqqCp1OR9++fVvtGeMtfK0/0P77dPr0aW6//XbOnj1LfHx8y3XKQkxMDO+8806Tn8etKnzdiW/0JAayqsZMzrGzjBvSn9qa6mYcYSAwBbgNGEx9xcBF4BRQApyzvDcAnYDOQAzQE+hn+bvLsl8p8BmQCWwAmqcOt9q5Pv74Y8xms8/N8n0xflX2qWVw9R12797NsGGuQ/W8jW0CDnfHg69df1/rD/h+n1xVHGoLc6PJZOKaa65p8n5XjMNVUVkFmbnf882p83SJirVsdScVYjfgIWAfsB9YCKQiVrAfA4sQquUEIBShdh4BjKFO7Zxu2RaPEMKDgGnAS8ARwIgQxP8FioCVlvZNZ/To0ZSWlkrHK4lPMnDgQK66ahyQAywFbkeMm5ahvLyctWvXag9gGRXQMXCW97ut/Hxqa2ubFe/b7oWvqqq88f52+iZ0Z86o3iy442rO/PCd9dMG9uwLrAK+B/6BEJhngX8i7LnhwK3A08A7wPFGjmflMvA/hKD9JdAbsRr+NbAbCAPmAJ8jHlBj3f6uIArLq6oqqx9JfJayskHAj4DHgfcRmqIvEBPZ/l4/38iRI7XXshjDlYV1MmVV6bqqOBQUFERGRgZ6vXOfG8XFdveJAkYhNKO/BBYDy4FAevbsyZ49e5p8RG94O08W/1Ine3qsplJeVcuKrGx+PvFWKi5fRO8fQLCxayN7DUOofvOAnyMcpN4F7gZiEULyA6DhogpN4wRCqA9HqKOfQTyQfmQ59/+Apl++rl27sm3bNjnbl3gFb43lI0cWMnDwH+h21Vp0+j2AGaFNWgJ8gxh7TwLeCQHp3bt3q5WBk7Qu1snUP/7xD8B5xaHk5GQA7r//fpc5HNQm5XbwA25ELLw+AAoRGstPgbcQz/I/AQ8C3fjuu+8YMGBAk7+bxzbf3NzcLCDL0+M0lVOlFSRG2w+s2uoqLpaedbFHH2AZMNHyvgJ4FfgrkO/yPF2i4qipqebC2TpjuqLTo5qb60B1CLEieAqYC8wHBgBrESvhhxEOWo1z9uxZhg8fTteuXSkuLuaRRx7BZDLZ2UFc2UYkEke8NZbDwuCvL01h4m2jMNeeBYKAm4FJwASE1un3iJXwB8DLwHs01ynRNkucbYWwKVOmyGIM7RRHFfLq1atZvXo1iqIwbdo0KisrAWF2OHr0KOCpo1gA4t68G+Hz4zhpKwO+RgjiM9pfSIiOoKDm2cd9Xu1su6qz5ud899Nd9QSvayKBfyEu3ETgEkII9wR+RUOCF+Dc6R+E4FUUfvPiOoaMnuCB4LXlEvB3IBGhxjABNwC5wL8Rao7GMZvNFBcXA7B+/XpycnJISkrSPpdqOElrYzQaueW6a7l4zjoRLgc2A7MREQOjgTVANcJ3YiMiOmAO4iHYfKwrn9raWt58881mhYBI2p68vDzGjx+vVQAzGAzaCnfHjh1MmzaN3r17e6EcZ2/gL8APiFXtVITgPQQ8h/DVSbBsuwFHtfOFCyfcqvHuDJ+vamQrPLZv38758+eZfMcoN/bUAb9AqA6MQA1ihv1HhAqhaQxJG0dsYl8CPEik7ZwqhGPWG8ATiJXvTGCc5fXrTT6iMwcAX6qSJLlyadzppRbItvw9hLjXf4mIkV8J/AGhjXqJ5kYFWImLi2PHjh0eHUPSNsTExHDkyBEtwqOiokJb4R4/fpzjx497eIbrEPeabS6I/YiFzzs4Lsr8AgzUVDm/H5sbR+yzK19nNhz3A/gHA7uAFxGC9z1ErO4vaI7gBdi7ZRML7riane9lutXe2C0aY7cY/AMNDBh5G+HRDZfkg/PAYwib8PsIh6//WPre2L6usTogyCpJktbAapNriE5hXSyvzgJ/Q5iE7gUOAHEIjdBhYAbuRSw4xydyuEuajPXZ3zIJW4YhtDC7EYL3EvAKMBRIQThR1deGuhK8nmRV81nh62hYd49AhDPTHsTFPImwM92JGMyth06nZ+F/PmHphr0c3PUpJUXfu7lnAeKmmImIJR6DcFK5t1n9sD6AysvLZZUkSYtjze3siiGjx5N49RCieiTZbDUjYuAHIjQ++xHOWK8jzDCjmtWXoqIihg8fTlFRUT2vWYnv4sypyvNkGd2BNxGLsjuAC4hQuJ4I35vcBvdWdM5FpSdZ1XxW+MbExBASEkJFRQWK4s7sdwjCUel3lvd/Q4Q1vO20dVDnMBSd91JBar0YcQO3TJhCcr9rCDX4A/D4vz8iZdRY/AKc2yfCusWweM1ndImKs9n6H0T/NyDiJP+LmKF1alJ/9Ho9t912GzNmzLCrkiSRtBQmk4l7ZvyM/iPS7LZH9UyisvwyMxcvp1v3BEaMvZffvLiOgCDbe/pdxFieiQgDHIzwMn0dEZPfNIqKioiPjychIYEdO3ZoXrMS36WuOIeo7VxeXu6BmjkQYc47hLDnXkaYIuMROR1cOejaYy0Va61CN2Dkrfx4+gwGDhzYzH75sPA9fb6Sfd+eILJHYiM6dT9E2MIu4GrERf4Rwov4osu9yi+V0a17T7ttQcGhdOnmnqOTK4Zc2593Ml9n2wcbmTmiB3NviGdm2kAiw4Koqap0uk9ZcSFLpo/i3OkfHD4pQjiJzUM4rcxGxEu6l01FURTMZjPdunXj+PHjLF++HJCJCCQtS2ZmJk889WfMtWZNwI4Yey/d4hKYuVjcgzMXL2fS/YuJTexL75Tr6RIVZ2OaMSMmn30Qji3lCBX0IWAWzVVFq6qqVUcKCwvz8FtKvIntM6mwsJD169czbdo0tm3b5sFvdT3ClPE0IvnRWwhP+0WIUE/30en0hHWL4eEX1jJi7L2Ya838/sklZGa6Z4Z0hs85XNWaVd7/4hBTRg12w5Adj3BUGoEYsH9DzGbcc9QwnRC6/Z8t+jufvvkS5RfLGD5sGNHR0Rw+fBiTyUSnTp0oLS2loKDAqQ2pc+fOjBo1iqCgIHJzc+utLgP99XTvEsRHm9a51SfnrEQk5chErIZ3Az9D3Ez2aI4Big5VNWMwGNi6dSs//PADS5cuZcWKFXZObCtWrPCgXxKJc4L89cz9/d/QBwRRVWtm0v2LXba1CuTXljxEnyE3MHzMFHZtXsuFc2f44dirnDv9BsJ/43YgAyGIfwp85+qQjTJlypRm7yvxLtaa66dPn9aiMs6dO6eFFzUdAyKm/BHE+vIb4AGEBqV5mM21lBUX8vzDP2bZpn3NPo4tPiV8Sy5Xs2XvD6z88zOoqkqXqDgnq0Gx9FfVexAekWEI2+4MYJvb5/ILCCRl5K38YcnT3DggmYBFv3TazmQysWTJEvLz89HpdJr3nTWecPr06Q0KsMa8PxWLkIyI7UHZWRPVLuugfo2wY/8ToZLLRKjnHkdMPASaY4Ba5yVoLall9Xi2Ij2gJS1Fz66dmDwggsjISGpqzZSWV3PEdIkjpouUlTvPu24VwoAmrF9b8hClxdmo5jsQasPlCBvwAUQ0wKvN6p/1wS7v/bbBmn/giy++0GJ2wT5mu3kMQZgo+iE8659GJMSocrlHWEQ00fG9CI+K0yZ+B3Z8RFCnYEqLi6iprkLR6Rh00x3cOXuBh/2rw6fUzqnX9uNnI5M1j2JngheCUNWVCBtoGLAO4ajhvuAFqKmqZO9nm0kf3JsAv4Yvg8lkYs6cOezatYv4+Hji4+PZuXMn8+bNa9SO2pjj2Jhx47h9yn3U1ta6VEvXcRkx438AETr1KMIzuksD+9gTHByshSFJD2hJa+Cn1xERHMiIxHB+MrwH9wyJo09t/xKWAAAgAElEQVRUCIob6uOZi5ez6PVs/AODEBPOqxF5QEKB/0P4RDTdFgzY3fuFhYXceOON3HjjjdIc0wosW7aMnJwcO8HrOQ8gNIT9ECaK6xEqZueCV9HpUBSF/sNuYvaSlzQzyKT7F/PH/26nz+Drqa2pJiAwEFS4Nj6aB8emMvtH8cwbmUDnAM98hnxK+D703L/pHBbeQIteCNvuLIQdaB4iLeM5QIT3LF7zGV1j3Etbd0t6ulvtMjMzWb58OQMGDODQoUMcOnSIAQMGsHz58kZ1/rbOA9aA8cTERCZPnkx8fDx+Cmx4/RWGDx2sGfOd0TnMVsC+gEhUYELkn94FJDnbrR4XL17UVuKVlZXSA1rS6kSFGri1fyQ/HtqdpG6dG2z7+PgUlkwfRXWlNd3rWUSigxmIamETEN7RNza5H2+++SYHDx4kPT2dxYsXs2fPHvbs2SMT0rQgtiGk3sO6CFuBSNLyAiJsqOF8y6rZjKqq7NpcZ77r0imAwT2MTB4cR4S+nLlz55KzfTtz586hrOQMIQY/ggL0jS7Y3MGnhO/O99/mUpkrQ/gUhDfzAETY0HUIW6hA7x9AZfllQOjnG2PGjBls3LjRwx67h+3Ked68eVx77bWsXr2aQ4cOacL77ay1HDt2VBPQjlwqO+ewZRsiX+5+RJaWXQhHM/cYPXo0UVFRfPdd8+1mHQ1PHNUKCwu5cVSaDHWxoWtwAGOuieaeId2J7+pcCA+88XYXe69B1NjehsjJno1Y5TTtkTZmzBhycnLsbItWU0xYWJh0TvQCttfQURNozUPg6rnXOAOBLxEhpWWIxdgDBAbpuar3taBFyoj/EbE98A8Umj//QAMpN4/lmawcUnt2YdrQq5gx7Cp+lNSVmDCD3aLLnYVWU/EJ4WudDX3+/lonn/ohgu7fAkIQqqehCBtoHbXVVZRfPM+S6TdzznTK5bkURaF///5cuHDBa/1vDHd/xJiYGH784x834cgngZGI8IwIYAuilGHjbNmyhcLCQj7++OMmnK9j40mqzl/+5lH27PqcZcuWtUDP2jdRoYGMGxDN1FSxElZQeHx8CgvuuJq9WzY1sOf3QBrCrqcgnGw+QIwFz5kyZYpMz+oFbK+hYxhRbW0t/fv3Z9euXcyYMaOJWrgpCDVzImJhNhixAobK8kuc/PYr0Jx2xf8zp05QU1WJX0AgNVWV9O8RxYNjBjMiMZyuwZ6lNm0qPiF88/LyXMx8ohEz2ocRevsHEMkmGhKcqs0FF+h0OpKTk9m8eTNz586lV69eXp/FeIuLFy/Sr18/FEXRZoWuArwteyDyj65AxLStQWTKcg+z2SyrvzSCJxVzrPt+sEFMLNetWyevtwsiQwIZc0009w6N46V3ckgZNVZbpSg6PT37pzjZqxax4r0dYYa5BZEwwVnbprF69WqnpevkKtg9XI2bVatWERkZyfr165k3bx5XXXUVjzzyCE899RQhISFuHFlBFKZ5C5H34FWE1q8uM9XDL6yzu3+sq9w+qSO5afw03nr3I+bOncvF0rNu5pHwPj4hfPv166d5EdcxAjGbGYlIen0TQpffdMxmM2lpaaSlpbWI+sCbZGZm0rt3b+bOnas5dV1/482N7FWLyJP7EMLzeRlCGLv/8zY3P2lHwFnGHXcc1QoLC106lMjr7ZqI4EBmjBrA1fFR2ipFNdfy3cGGQjw+RgjcXYisRTsQdmH3CQgIIDY2VnsY6/V6bQJsm9hfroLdw9W4mTZtGkVFRUyYMIFZs2bRs2dPcnJyiI+P58iRI40ctTMicdJChNPpQ4jSsPbj7B/3383+re/brXKDg0N5PTOLzf/NYMyNw9pcFvhEqFFeXh4JCQk2W36BEB7+wFZEiIFn2ZkyMjLaTUyr7Q1hTYyRlJTEDz848/62ZQVwCliN0BJEA/fheGM6Y+DAgRQVFUnnKyc4qsrcdVRLTk52KWSbmw+2I3Hh3Fnmzp3L/73qOpzI0DmEissXLdquU9RN0ucgQk6GAL/FnXKFVVVVmp0X7HND2yb2lyF67mEdN9YsheXl5bz55pva59XV1Vx33XVNOSLCxDYYkSTjHoSpzR5rWNDlC+cJj4rjhjun8u22DZgvnWNAnO8kV/GJlW9MTAzTp09HqE0zECUA/RG23nQ8FbwgVr/tWdWXmppKr169xJsG1SRZiHqUZQibyPsIW7lrgjp1Jjc3l0ceeYTIyEgOHDjglT5fSVid5rZt28acOXMaDDGzqttcJfaPjY2Vkxw3sPpKfHv4MPfcMxWd3j60Y8jo8Qy66Q4UsEndWoXI1TvX8vphRJUa1/mmbcnLy2tUKyFD9NzHZDLRr18/AIcFVlO5BqHVGAwcBYbjVPAqOlBVDJ2Cmb3kJX73p2f53bRbyHz1Zda+VT8pUVviE8IXYM2azxCei7MQ8azTERlKajw+tl6vb/eDJTMzk+pqkZwgJk6k4dPpXSkutiLU9acQRcw/pSEnlPLLlzCbzaxfv57z589z3XXX0a1bNymEbbB1mnvssce0uFBn9r+te/YTEec63C0goHUdO9o7MTExhIWFYq6tRWfJx27NE32xtAQUxUmM/CsIZ6xiRCL9nQjHHM+pqKjgs88+88qxrmSMRiMbN27k4MGDqKpKQUFBM480GshBFNv4HGGSrFNP+wcaCIuIYkj6BB5+YS3Dx0zlclkJY66J5varowgO9AkFbz18QviGhSUgnCSuA44jjOdveOXY1kxU7Tme1bqSsiYXL/xeZKwy11omJoqiPZTq+ApxHY8gVG85wFVun/PChQt2KiEZclHHsmXLXMaEhoWFMWxAX878cMLl/sePH5dOV03EZDIxb948du0SfhBX9+3LL/70gpaEw9a5po4diGfK19SlZb3B476oqkpRUZG0/TZCXl4eSUnu5R9wzb2IEoBhCAer0aCcJSxCPMv9/AOE1/KwUdw7fymxiX15fMmf2fr+hkZjyNuaNhe+RqORyspChMfaJ9TFrrpHVFSU5hiRmppKly5diIiIaHImKl/G6rhgzUxlxWAwMDx9PEPSxqGqZvwCAh08944jHjb7EUnqdyASi7uPwWDAYDDw29/+tsOHXDhLEGAbE/r9uXIef1VUsHLOrViHXHvXxLQ2juF6mzeuY/p1V9E7KpjQ8G4EdursIkPccUSmo80I7c8nCB8Sz2mK13tHw2g0kpCQwLFjxzw4ykOITIYBiLz99wIVoKpUXDoPQP/hoxg+ZioXzp1Fr1MY3bcbt/aPxODv/Yp13sZj4Zuamjo5NTU1PTU19dHm7F9XfNsaLuBeiSd/f1GuLygoiEuXLnHp0iVycnIoLCzk+++/b3ImKl/G6rhQVVWleV/q9Xqqqqq4umckutpKho+ZqmVssceEyIVrXfluozlhGOvWrWtymM2VRl5eHuPHj9d+A1vunDCJ974qonOXCAI7OZtxzwU+RDgBKT6pifF0LLc2QQF6busfxU29I7hUWsLwMVOZszSDiFhHlf8FRJ3gFxB+JW9SV3rUg/NL2289rBqybdu2MXXqVC2EtOlJNJYB1vKPv0VUqat7tlkTKh3Y/hE733uTw7nbmTw4jv4x7tn2fQGPhG9qaupkgNzc3E+A0tTUVPfyNdrQr18/i2A0445HIkD//v3ZsWMH8+bN86ieYnvC6vAzcuRIevfuzciRI5kzZw5nzhSzdfMG/rDsrzzx2scu1G9liFXX+4g8uJ/SlGxYjlgfOB1BFe34Hffs2ePUkWrdW//loVv78di4QRT/8B06vZ9NStCpCCdC6NLlG+bNm+tzmhhvjOW2YkBcGO+8ncVP5j9J75QRJA8a7qSVGREB8Ijl9TPAy0DzV0jl5eX4+fm5nER1hPHhiDWhRkZGBqGhQhAaDAbMZnOdw2iD6BFOt48h/H1+AjxXr5U194F/oIEbbp/A1wcPEhnivF66r+KpJXooIuUUiAjnwQi9jlOcpdbbvn07Y8eOdSvt3tChQ+nbty/FxcVER0ezcOFCl8f1Fr4SSvD8889rr0tLS+1WnmfOFNPfCGqvaN63xLQpiuKwCi5H5MFdjXDR/wiRku3DJvVDURT8/PzQ6XQ88cQT7Nixg0WLFvHoo763WPLGb7d48WLtO4Iozu7n50dNTQ1BQUGUl4ucw4pOh2o2U1tTzfFvvsRcW0NVRQXC2ed1QIeiLMLP7xUWLtwLtOx92ww8HsuOtObY0QFpPQO5ut+11FS5rmAjIii+Q4yDuUAcYjxcbtZ5t27dyldffcX999/Piy++SGRkpPaZ7b3jylzjK88XW5rTp169ejmtTqTT6diwYQNjx451I4Y3AOHrczfi95iMWDDURzWbNXtvfEQIQQH+rT6ePP3tPBW+jrrHrg01tr0xbbedPeueqvmrr75i+/btbnfOWzjrd1vj/FqCUVfJ9WOnct0dU9jy35Uc2PGRTcavauDHCDXcLGATIh2le7WGe/TowR133MErr7xilwvXl8uzNfe3cywFaft9a2qEo5tV8IJ4GFiprRFe6dWVQxDX1h/4C6r6NMXF0Lt3b5+7TnhhLHvSzlscyjvErPt/w+fZH1JdWYGi06PWy/W+HuEJ/Q4wFvjM8r+4yec7efIkQ4cOBYTAWbFihdN7p6Hx0V6eLw1x6NAhHnvsMTZt2kR5eTmBgYGYzWZeffVVnnrqKV7+z1oe/MXPuXzB1X0fjEiekY4omDEW4dnsiIJ/YCAhxq788/9e59NNmRQVFbXZNfTkvJ7afEuBhsoQuUV6errLknu2VFRUdEhbY1N4d8M6/r3yn/Ts1Z/OYcZ6qTaFym02woEhALHY+albxz5x4gTffPMNu3btsvNiDAoK4q677rqibF/OEsA3zW41BJEQIAgR9iI0A7Gxsb56nbwyltua7nGx9O4eqWU2QjUzZPR4rhkx2qHlLoQjVj5i0f857lYGc4XVH8JsNtOtW7cOV7ozOzub8vJyFEWhsrKS6upq7r//fnJydvDLn09rQPB2QShZ0oEiRKIUe8Gr6PSk3DyWxWs+5e/v7Sf3f99w1+jr27U/j6fC9wvqZsyJiBxvTWbr1q12qwhXdIQb2BtcmxjLw7f20+oiO2c+8HuEjeVV4EG3jp2Tk8N1111n58VYXl7Opk2bUFX1irFxOUsAXz8FqiuuRaj1reERv9A+OXXqFPHx8b44ifTKWPYFTCYTc+fO5Y2NH3L92Kl8mf0uX++sn5BBhOFdj0hjm4yIBR7a7PMaDAbuvfde7rnnHoqLi6moqGhSRrT2bCNetmwZxcVCc2Br7hIrfZXqygoXe0Yj8hIMAwoQ0Rkiv4A1j4Ew6dRi6BRMYo/u3DMkjlijo19L+8Mj4Zubm5sFJFqdMyzOGk0mLy+PuLi4Rtv5ooeoL6KFJjWqTViCcOcHWI4Qxk0jMDCQ5ORkbrrppiuuAoxtVitr6FrK0GEYOgU3sFcfxCw+HKHWn47QNtjja7mdvTWWfQFrWNL4m4fz31df5pm12wiLiLJpoaDo9Jb62adBSUP4PlidEV2VMWyYiooK3nzzTTsTRUVFBWaz2a3Snb40fmwnAg1NCjyrzxuPiMK4FjiIELzHQFFYvOYz+g8bxYix9/Lw82tJve1uzJdLuTslltAgf0++ms/gceqP3NzcP3t6jJiYGMaMGUNGRgYBAQGa4d7qNBQbG8u4cePa5YywLdBCk1wk9bdnBcIbehXwJEJo/AZbt35HdHo95tpa7bcqKCjg6NGjfPrpp8CVk/vWVp116NAhvj19kcG94lzEk4JYMG4BIhEP83twlqEtPj7eJzMkeWMs+xqRIYHMuy2Fj29MZ8vbb6D386e2pprI7vGYTuZbku6fB+5EeNnORNiCZwOveXTuoKAg4uLiOHbsGD179nTZztFG7Avjx3EiYH3tmB8/Ly+Pn/3sZ824n/silCvdEQmWbMJMVZUl00fh5x/Ask2imMaiPyzhtpSENqtA1BK0eZINKyaTiR9Pm8727du1Vcbu3buZN28eQ4cObde6/bbAumrbvHkziY1mmXkNkQe6ErES/g8NzcvMllAba/mvmJgYJkyY0OSqP76I4yzf+n5D9k7uGns7P/vDi+j8nF2bRMSqKc7yfyINFbSQGpzWw+Cvp3PtRcZP/Qmzl/2fSD948TzDx0zlgb+/QVhENAGB/gjfh6WIe//fwBNNPldSUpImIMrLyzl69CiqqtrFxzveY82tmtUSGI1GevToYVcG0LEkYFhYmNb/mJgY9u93PymSYAiwHSF4tyKc3+o73aqqiqIopPXpxuDuwVeU4AUfqWoEYpVhMpmIjIy0u+msVX0kTcN2ojI6LY38RjPNvA2MATYgSrF1QQhk17Z4q5f6999/z/fffw8INXR5eTlVVVXMmDGD1atXtytBYzvjX7FiBcuWLSMnZwe7J95BTVUlm15ehrnGcTWbhBC4VyHUaONxdd2CgoI6TGy6L2EdD58eOE7y1YOZdP9i7bNFr29h3fNPsmvzW4hSdacQGqGnEQLifpyZDpxx7Ngx9Ho977zzDr/61a+0lLBBQUFMmDCBZ555hqVLl2r32KJFi5pdNaslyMvL4+GHH+ajjz6ivLxcSyhTW1uLoiiMHz+ekJAQ3njjDeLj45txhpsQmoUQ4D1EOFF9e3DX2B48+NxqxlwTRWJEZ0wmVzbj9ovPrHwlLYc1L+47775HzFXxDbTMRsxCzyBc/a22S/eYMWMGd955JwAbNmxgx44dbWq/On36NDfeeKPLAgi2uCr8LWxZwmFEVVVOn3CcxCQjZu/W7GF3ABddnqe8vJyNGzf6osNVh+Dq6E6k9+tWbxV1zlRISHg3grtEAC9SJxR+iQgXazwaw4per2fp0qVaFj69Xk9lZSVr164lPj7e7h7r0aMHRqOxSVWzWpKYmBhCQkLsHA1ra2vR6XSoqsrGjRtZvXp1E5wPbbkT+AAheN8A7sKZ4AVQzbXcd/MAEiN8Oz+zJ0jh2wGwOqDckj6asbeloyiKliGmPrmIikgnEJ6g1moijbN69WrWrRMxw2ZLqsuG0lFaKwO5Ixybw/Lly10WQHDEUfUXWC9LmDP6IARvnOX/GBoSvLb4msNVR6JvdAh3XhuNn80YCI+K5eK5MyKBjN6PupjTEoSQ+IRGQp81qqqqyMnJ0ZJKWD3la2trGT9+vBaCZDAYuOuuu9i2bRtnz57l8ccf94l0uGfOnNEmAlaaJ2xtmYG4pgbgn5b39hqkqJ7JJA0cxo3jpzFsyOArwqO5IaTw7WAUFxczZ84crh0ynIi4eBetDiEE71dAP0QIRvNVpTqdzqn9qqHqQJ5gXcXaep26mgRY7W+KomiqP/+AQKqqKomIc+0kI/JjbwdiESrnMcAlt/qXnJzM4cOHm/itJN4kvmsnJqXE8PiEwSy442p2vpeJqqqUFhfVVQtjByIN63eI8bCT5sYCJycnA/Dxxx9rzlUVFRUEBweTkZHhdS9nT8KWVq5cyWuvvdbEQvcN8Qgiy5sfQpX/a2wdOgM7BdMndSS/fWkjv1v+Omtfe5m316310rl9Fyl8OxiZmZmsWLGC3du2sOb9HAanjXPR8gfECvhThIDZBtzSrHNOmzbNzn4VFhbmsjqQN9Sx1gIItkkx9Ho9EyZMqDcJsLXxnvihkJHj7uX+v7/B4LTxlJz+wcUZRiKyInVDpL8bS1PSE9bU1LQrO/iVSlSogS//9zXD0l2NARAT0eHAl0AvRHIOZ7mjG+bo0aMA9fIZrF69up6pwxtjwNOwpby8PC03s2c8S11u5ocRBXTsqbx8kdlLXqJr50DuTokl7AoJJWoMKXw7KIqicFPvCDrraojq4Wo2X4YIAXgTCEWUZZvXpPP079+fCxcu2G2bMmVKvXauhGNziImJISoqyk5VVltbS0hIiCb0nNl4P9z8Lls3vsG6558EVCeOVSBWuB8irsebiHzZ9Z2r6ooqCEJDQ9m8eTP33XefdLjyIXonXEX/niJFoIj7FfgHGki5eSwPv7AORWcCbkQ4CEUgfCMmeb0viqJ4PAZc+S64K9ALCwvp0aMHCQkJnD9/vtn9EKvcVxGZ3aoR8e71nWfDo+LokzqS6FADd6fE0NlHC9+3BFL4dnA+fvdt+vbpQ3hUHKFdI+kSGevQogqR/9kagvESIjVlw7dOr0EjGJp+Fyd/OEVhYSFFRUVO1cFWamtriYqK8tqK0GQyERsbS9euXenRQ9isc3JytM8dV8f+gQbCo0SilxOHDrB3yyYnR52F8AYPQlTEmY54sNTnUtk5AAICAgAIDw8nLS2Np59+WobM+RglZ84wZ85chqeJ5BqKoqOmqhJDp2Be+M2PLXm7LyEmWi8hfv91iMo73sPq0ORJJixnYUsTJkxgwIABjSbMALFiBggObiiRTGOEIFKr/hRx3cYhHKzq0yf1Bv7w/GtMGBhDYDuowetNpPCV8OE769myaz+/X/MZT7z2sRNVtIoIwfgpQhj/BtiIWP0558j+nfgFBnKhrJQ9e/bwu8VP8s03By21m53TvCw5zsnMzOT06dOcPXuWEydOAHD8+HFtFVBrMLL3qzxtdVxdWdGAmllBTD4yEEUSliJSRjbshDJv3jxycnI6VOnL9khmZibPP7+CyGB/xt7zEx5+Ya1WoP3xf39kkx2rFuH9vADx2y9DxAMHeLU/ruKBrTSkUnYWtvTtt9+Sm5vL0qVLXe7rmKnq4kX3HAfrE4fwhbgNUUv8Zhwrpw0YeRsDRt5GeFQcNRdLGTcgmgC/jieKOs4aX9Ig/WJCMPjr+PCgiaqKcsKj4ygpchRGryHyr65HhA3sQSSTyHN6TNvc0pmvv0rm66822Ie77767+V/AQmFhIYmJiS69iW8bM44tH39IckwXp5/XJxDxvacivDN/iRDCrtHpdLz77rukpaUBMla9vWDVSOzMLyE2sa+2vd91N7Fr81ta2Ui/gOfxCzhNxcV/ITJiJSHGwRmv9OO2227j9OnTDBs2jOLiYpYuXcrjjz9e7762zYSVl5fHfffdx+rVq7WwpVdffRWz2UxeXp7W3nFfEDWq8/Ly7KoSNY9rEaap7sBhRNhdgfZp19ieRMT24L4n/gZAz66dGHN1FH76jid4Qa58JTYkRHRm4qBYfrnkRWIT+xEUHOKk1TZE8vn/IUJt9iBiIhtH7yccKQJd5EZet26dW/Yp2xClAwcO2K0OrGqzpKQkLaTDli2ffMT9f/8vKaPG4m8JJ/ILcFWEOw7hWDUVOI+w9zYseEGEZWzYsKHRdhLfZERiOCOTI1AQscAXS0u0HMNRPZKpqaokLPxzhOPdSURO4lyEB7znfPjhh+zfvx+TyaSF6yUkJKCqqp062NZGbLuitYYWHj582E4F7aoq18yZM4mJiWHt2rUeCN5xCO/w7oiV7/XYCl6A3inDmb3kJQCSI4MZe010hxW8IIWvxIGo0EAmD47loaX/xC/AQFTPZCY98AchoDSHlAJgBLAGUYdzLfBXhErWNdY6t5WXXau0unTpYhdf6AzbEKWZM2fy+eefa8J25cqVqKrKsWPH7PLlWqmpquQf99/N/q3vU1NVid5SkNu/XlzvKISH63BEqMmPaEqhH296rkpan0FXhXH71ZHodQozFy/ni4/W8/df383pE8JrWSRb2QdcR6eQPKAnQvjc16L9slUHW23Ejok7DAYDQUFBKIqCn5+fXX1dZ+Tl5WmlEJvH7xC+ECHAfxFRESWAmGjPWZrBiLH3cuGcyIg3oHsYt/cX17YjI4WvpB5hQf5MGRJHxvt7+O1LGynMP0RtdRWBdivJckSg/IMIdex8RBxkL4/Ofe7cOZYvX056ejoHDhywW+EGBQXVC1HKy8vTEhjYouh0+AUEoOicO3GoqkgCYi20bl/y7LeIpAqRCIE7BPjaZZ8VnY47x41n/PjxPpGfV+IdkiODuWtgLIF+eh7/90d22hL/QANdY3ugKKe59kdPAysRjlj/QXj1tl24jDWz1tKlS9mxYweAVqzGuxiA1cAzCFHyBMI5s+5cnUPC6J0ygkn3L2bm4uWMSOzKTb0irrg8zc1B2nwlTjH465l5Uz+71WNlubNY1ucRquc3EEJqH/AAIsygeaxZswbALsh/5syZgPDCdMcZRDWbqamqAhoe5GY7oR2FqO401vJ+KbCYhhyr/P0DqK2tITYmGlVVfSI/r8R7xBoNTB4cyzsHdAR26qxVtKqurODsKeHIt/uDNQgt0F5ETugHEfVp7wWOt3qfq6uFhsmbDoz1SUJovFKAC4iJeF2EwJylGXy94xMunBN2cEVRGN2nG/1inJmyOiZy5Stxift1gXcjBuEaoDPwf4jk6Vd5tS+qqjoXvA3Oot1N43g3YnU7FjiHSCm4kEYT6itouXh9JT+vxLuEdw5gypA4ai6VMnzMVPoNG2X3uTUmePGaafz8yW/wDzQhhO8+WiIeuClYV5iu7L3NYxJiopECHEXYd+sE75DR4+1Wu346HWOuiZKC1wG58pW4xBq2UOGWE8Z5xOz3A8Rq+E5EBZOFiET1nuaGdU2nECOXz59r2k6KAqqKSJrwN+psdR8CP0dUtmmYzp0788033zhd3UoP5yuLTgF6vt69zakfQXVlBYZOwYSGd2PZUynUVHdGaFAmIuKB/4UwZbifBc1bWL2jPc/NDCKk6llEpiqALETsu30yji+z3+He3wrHx0A/PWOvjSLO6H5hio6CXPlKGsRkMjFjxgxuu+02l7Nn6+y6szEcYQPqhxiYIQg13E7E7LhlaLLgBVB1iByz3yIE72XL+9txR/ACXLp0ifj4eIKCglqkMITEt7BqggINQpAoOj19ho5kSPoEzZlIcA6xOnwAYf/8JXAA4R3tPoreO0knvGNfHYTw6H4YEev/IKLkaJ3g1en0hEVEs2j1pwB0CvBjUkqMFLwukMJX0iCZmZlkZGRoWaIcCTAEEdUzGb+AQC6Vlli2FiEG5rd4SkwAABi8SURBVATge+A6hCfoekR4UltzM0Jt9gKibvGHiMIR/3TSVhGhUY08wNqydKKkdbBqgqqrKgkMDATVTHhkLPfOX8rMxULT8fi/P6JrrHWsvIC49/cjbKSfAX/H3fKEqoMTYXPxrIKWHqG92oOI4/0WMYl43sl5zPQfdhOh4d0IC/Jn8uBYIoJdhfFJpPCVuIXJZHI6g66qKKfo+BHNEcWeTYhV8JOINHMTEXbVVcDVLdhbV9yKiFPORgjbAoRt93aE7coeRVEYMfYenlq3m2ff2sb4SZO14uK2NFY6UXLlYLXrb9++nRk//TmXy0rsPg8N7+bgxHcAIYCfRJheHgYOIsaCr3MdopDEUwjv7ecRdt49dq2s4UTWrGARwR2rQEJzkcJX4haZmZkcPXrULmjfPS4Cf0QUnX/Jsu3nCCH8IULwtWTYQRAiB/Nuy/lGImIQFwL9EWkynaOqKrs2v0WP8E7MvW0wUV27UFtbqwng0NBQ7bUMLeoYWBNYDBgwgFf+9SIfv7Oe8M726SXjkvoxYuy9WlIZkf/7jwhhtg+IR2iBPgL60lT0/t5NZ1mfbohkMruBVERt79Ho/R7Bmd3aNpzoiedWcndKbIcqkNBcPBa+qamp6ampqe5nH5C0W2zzxgYGNlWdVISwffVDqOMuIVai7yOyBP0DYRf2hiD2Q1ShecVy3tWIB58JkRCgJyKMqL7zjC19h45k9Se5jB8QTVCAHpPJxLx589i5cyfz5s0jPDwcVVWvqNAiOZ6bRmiQP1MGx9E7qi7z1MRfL6Lou6MkDxouKoZpGqN9CGH2K8QE8BbEyvhlwLlZxxm11VXe6r4DnRF5q79FOFJVodP/BbgGyKbWaZUvOGcSPhL9YkIY30HzNDcHj6cnubm5n6Smpv7OnbYmk6nBz0tLSz3tjteRfbLn5MmTTJ8+nenTp/PAAw9w5MiRJh7hKMIR5ffAHIRAjgcesvwVIuzDOy1/XyFWzw0RjbCpDQNGIwSvbQrLnYj8zK/TFI/Tnt1jGdmrG8XFxQA8/3ydnWvhwoV89913jBw5kunTp7NmzRpOnDjRLu9xW9wdz419T/DN79pSfUqJUOik6tn93QXef+15jn+zlyG3TELR6bmq3yBSb5nIS7+9D6F6/hfwFqKw/GxgLqJoyauISeGJFumja4IRzoa/RXj/g5gUP4S5VoxvRVFEcIBN6J6i6LjmhlsZ89OHSIn0o3c4lJz1Tn5rR67Ee6lVdQORkZFeadPayD7VYZuzuH///lRXVzNkyBAul1ewe+8+qquruFhaYinD1hDngD9b/oYC9yCctHoickXb5os+j/BALkQ8vPwQNqgQIBExY3fkEPA2QugebvL3DDQEYTBXNHidba/FqFGj3D62L95PTcXd7+CL37Wl+tS7t9EuFCn3o3UA+PkHMPXhP9En9QaO7PvcYhM+i6iM9XdEgfkfI2plz0aU43sJoZZuuRA9oYWahRD8XS3bPkfYpz/SWvkFBNJr0HDy9mxFUXSoqhlFpwNVJbxLF+4ffz3dQlresepKu5ekYl7SbGzr0ppMJiIiunHPT+fy7luvN/FIX1j+FiBsYMMRuaNHIGzFoZY/V/axs8AxxCo5G/gUIaibR7fIaPr0TpaxupIm4VgZyD/QwDXXj+bO2Qt4fHwKNU7VxYeB+zBGvkpZ8SxU1RolMAGRHWsNQhjvwTuCuDui2tBMRL5yKzsQdulP6u1RU1VJ3p6tAISERxAUHEqwsStJyb0IqDrfKoL3SqRR4ZuamjrXyeb83Nzc+r+SpEOj0ylszlrj4VEOWf7+bbPNCMQi1MsKwoGlBpFfugDwrkqq2FREsamIpKQkLl265NVjtzVyPLcczmrpWpNv3P/3//LKwjlcPl+KqppBUVAUhZiEvhQWHKLUlI2YOP4G+BlCFZ2IcAxcCBQjEtjsRlQUO4Bjcov6+CMmr/0RE9rbEfZbKxcQaWEzgFxNtaz3D3BpVz5/1sT5sybOFZ3k612fyRzNHtCo8M3NzW3JBKGSK4xjx47Vqwsa1i0aBYXyi+eprLhMUOcQyi+eJyK2B6VnTrsIU7Kl1PJ3sKW7b0dtbS0GgwGDweCTNqfmIMdzy2INRZo1axarVq3iyPHv8dPp2P3+W1yyhCX5BQRSU1WJqqrEJvVj1p/+yb8W/IQzp04gnAKfRZhjRiFWwOMQgvg+7Ksm/YDQ+pQgzDggTDHBiPj1ROoXeLgAbEF4+a9FOD4KHn4hi12b13KuuIhOwSF8mf1Ove/nHxDIrWPG8eI//ioFr4d4rHZOTU2dLP6lTs7Nzc3yQp8k7ZiYmBiysrLsqgyVFdtnfyq/KGbsZ061tmNJ0wgKCmLChAk888wzbd2VVkOOZ8+wNcU89thjJCYmkv3he3ZtbCebuR+t02zD9qgI88mniNjgvkA6ItPUQMQKNs7y1xD5wDcIk8xHCJtutc3nCtdcP5qJv15EaHg3Jt2/mEuXLpH1tyeI6pHE6ZP5ljSsIoNVTU01cZHh7d6r3xfwhrdzFiKXoEQCQHp6OseOHePkyZNUVlai0+mIjetOUVGRC7uXezSkDmsJrpTwoaYgx7P3WLZM5DdOSkri1KlTlJeX4xcQSIixKxdKz2r1pPV+flQ5rRhmi9UcY0WPELxdgHDLfxWxsr1o+V+Ao3e/cJiy3aIS0qUroeHd7NrNXLyc15Y8ROK1Qzn7w3FqLp+jR2wMffr0kalUvYR0uJJ4nY0bN/LAAw+watUqDAYDVVVVhAR35lRtDTqdTiR51wob2BMe3Z3amhrKztQN8ICgTkT3TOb8WROlxU0f+IFBnaksd992m5yczIoVK9i4caN80EiajNFo7/V87Ngx7bW5pprQkGBKi4VDYG11VTMnlLWIkCShPQrrFk3PvgP5Lu8rzpec0epUA/ToM4DKistcvlBGz74D+eFYHlf1Frbf77/92iEvdR0zFy8nqVtnbkyOINggRYW3kVdU0iI42r42bdrEnDlzOHz4MEWnT9M5LJxTRacpPH4URadDNZuJ6plEt7gEvv7c3venqvwyJw4daHIfwmOuoltcPLOXvMR//vhLSk9/z+lT31NV5fxhFxgYSHV1NWlpadqfRNJUHL2edTodiYmJ2oRu06ZN/HzWbPqPGM3SRQs4W3gSELG0iqKjU5hIUXr5fKlDqkrXlBUXcaC4yC4UyBrud7H0LI//+6NGjmBPiMGfm3p1JSHCWRifxBtI4StpEWxtX8uXL3catjN16lQ63XQTI+6YwltrXqP0bDEzFy8nY/EvOHvqO49twn0GDWPJX5bTKzKYB7Z+pK3G9Xq9nU0aIDg4mOzsbFatWiVXu83gzJkzlJWV2SXxN5vNnD/fmEdu69Jaffr1r3/NzJkzLR7EKiEhIXTp0oUHHniABx54QGuX+uZqLl26pLULDOpMp5AwAC5fKKPi8kWysrJ49dVX0en9MNfWoPfzp7amGkWvR6/TU1NdhaLTMeimO7h84TzffrnDLs6+5PQPLLjjavz8A1i2aV+D/Q7Q60iK6cwtg7rjr5eZqloSKXwlbYatgJ41/iZOlFzm5Lly/rJyDc8/9Tgfr1+jzeQd0en1LlcFabePIy4mEtOp77khuau23boaP3z4MHv37iUkJITrr7+e3NxcBg4cyIABA2RsbzMpKysjPj7ervBEdXU1/v6+lVy/tfp07NgxoqOjiYiI4MyZM1RXV5OUlOS0bUxMDOFdIygymaiqqqZrzFUAnC08iS6iGz+f9ysOmcr5eucWrhkxmuH/3979xrZR3nEA/+ZfUygkl2vTAlrV3jUSQzA1xg+TOk1Is8M0RZNYPANSpGqIlOZFWV4Qmb5IEG8aqrbyxhhoIjgdm1BUqZb7BoVKrb2tL8qknbFfTOIN2Fk1RKHgXAJam5gke1Hf1X/i2IntOzv9fqQK+/75h32//O6eu3ue/mdwJfRXxK98iOXv02jd1o7l9BK233sfBl85jYXUjay7p5HzvHEx929vw8EfdODRBzugp75m4bUAiy/VhZbmJii7dpjNXFPL32F4eBg3btxAKBQqWD6/8D78w0fw17+8h7Nnz+L69et4909vF3SBmF3sqbpWV1fXHPHpbpVdaIsNxwkA+/btMw8GepT9WF5Zxc30Mm4uLZtF+Oa1WXheehWel14119u2fTuW00vYs+8ABl85jX/OnMe3c1+v2ZlHevEW4n//EIOvnC74/Ac6tsOxtxPqrh1obuajQ1Zi8aW6ZBTK5557DsPDw3j66acxMjKCZDJpNhl3dHTg3Llz5o1RBw8e5JkrNbSW5ibc196KHdtacCu9gsXvV3KGks4vrl/+5zP8/tivzSblhdQNfPDuGcSvfIjVlRW0tm1DZ/eD2PXQnQOAJjRh/6574djbyYHubcTiS3Ut+2w1u/ACwMLCAvr7+7dUJxhb3RdffIHDhw/j/fffv6se4dqopqYm3LOtBfdsa8F97a34Ve9D+O/cTfzu/D/w5z+cwL+vhpFevFXQpNwhd6P93h3A6qrZHP2o+AlGXj2JPfe344GOduzp2M6Rh+oAfwFqGH19fejp6TGHM2xpaeEYug3m5MmTuHr1Kl5//fVNbyMWiyEUCiEcDiMQCCCRSAAAQqEQxsbGqhVqUbquo7+/f824HnnkEYTDYYTDYfj9/qodFO7tugeHVBlHf+GAQ30Q3y8tYVv77Z6y9uzqwmM9+6Ds2oFHH+pA29K38Az+BqGZS3hh6Ai6mv6HX/7oATyxvwt75XtZeOsEz3ypYaz1/PDd1glGo8p/9nVychKTk5MbbrXQdR2BQABvv/22OW1wcBDT09NwuVyIRqNVjXstkiRBUZSC6Q6HA4qiwO12AwDcbjf6+/sxMzNTsKyu65iamsLIyMiGP/+br2/g6NE7j/Fdv34dTx980Jzv+uCC+frnP/3xhrdP1mDxpYaS//wwHwtqDPnPvm62685gMAiHw5EzraurC7FYDIqiIBaLIRaLIRKJYGhoCNFoFLIsIxKJYGBgAJcuXUJ3dzd6e3sRj8cRDAahKAoefvhhXLx4EdPT0zh27BhGR0cBAJFIBLIso7e3F7Ism8snk8my4pUkCbquI5VKIRKJYH5+HgMDA0gmk4hGo4jH49i5c2fOPFVV191m/mN81JhYfKmh8A9PY1prxJ9atFo4HA6zOE9NTQEAUqkUvF4vxsfH8fzzz6O7uxtTU1Pw+XwIh8OYmJgAAFy8eBEA4PF4oKoqBgcHMTQ0BFmWMTU1BV3XMTo6ClVV17wDv5hUKgVVVSHLMqLRKC5cuIChoSFEIhH09vaira0tZ55R+GlrY+M/EVnCaLW4cuUKXnzxRXz55Zcb3obX60UkEsmZlkwmC86GDQMDAxgaGsKZM2ewtLRkFmej4EqSZC7rcrng9/vhdDrNaU6nM2d5WZYBAHNzcyiHrutQVRXhcNg8GMiPvdg82tp45ktElqhGq4UkSfD5fAgEAmbz71tvvZWzjNHsPDo6Cr/fD5fLBa/Xi97eXoRCIfT09ECWZaRSKSSTSei6DkmS4PV6MTY2ZhbkEydOmM3MsixjdHQUwWAQTqcTyWQSsVgsp+gnEgmzmAJAPB43Y0smk1AUBfPz8+ZyADA7O4tr167lzMvfLm1NLL5E1FCym5azSZJknqEa840mXOP9Cy+8kNPDlXFzlLF+9o1cqqoWXH813n/00UcFn6+qKj755JM1t33kyBHztcfjMWPK73HLmEdbH5udiYiILMbiS0REZDEWXyIiIoux+BIREVmMxZeIiMhiFd3tLISQAPRl3j6hadrxykMiIjswn4msU+mZ77MAZE3TggAghDhaeUhEZJO6z+dEIoFDhw4hFoutO42o3lV05qtp2mTWWxXAO+stnz+4eb56HBaOMZWHMZWnHmMybCSf83N5ZWUF6XQa999/X9bU9ori+fbb7wqm7d27Fy+//DImJyfx5ptvAgC++eYb7Nu3D4899hjS6fS628wekrIelBPPyspKyb+d1VSP++hWjKkqnWwIIVQAKU3TEustt3v37pLbKmcZqzGm8jCm8tRjTNnKyef8/4eFhYWcziKqodj2ZFlGc3OzOb+1tdV8n0gkcgZDAJAzaMGnn36KN954AxMTEwgGg2anHHYq9b01Nzdbvs/U4z661WIqWXyLND0lNE27nPXeq2na8KajICJL1Dqfb91aNF/n995UTR6PB6FQCIqi5PR2NT4+njMYwsTERM6gBSMjI3jvvffgcDgQCARqEhtROUoW37ymqAJCCK+maaczr/vykpiI6shWyWe3243BwUH4fL6CeU6nE5IkweFwIBwOIx6Pm2fBwO0hCInsVtENV0KIPgCnhBBRIUTtR7EmoppphHyOxWLm8H5OpxOdnZ2IxWLmQAfGYAjhcNicnj1owblz58zBFIx1iOxQ6Q1XlwEcqFIsRGSjRshnh8OB6elpAMgZ9zZ7oIPswRCym6Q9Hg/S6TQOHz4MAJiZmal1uERFsZMNIiIii7H4EhERWYzFl4iIyGIsvkRERBZj8SUiIrJYVXq4IiIy/PFvnwEAlpdX0NKyueP73/6srm+6JqoYz3yJqKHEYjEEAgGEw2GEQiGMjY2Vve7Zs2cRi8UK1tF1Hf39/dUOlagonvkSUcPQdR1nzpwxn/UFgGi0vP5AjI41hoeHoShKzjxJkgqmEdUSiy8RNYxgMAiXy5UzzefzIZFIIB6Po7OzE/Pz8+js7ITf788ZQCEajWJ2dhbhcBh+vx8zMzPQdR3BYBCKoiCZTAJAweAMyWSyYFsA4Pf7oSgKOjs7oShKzjrZHX0QrYXNzkTU0CRJwvj4OFwuF9xuN4LBINxut9m/szH0m9PpxP79++F2u82z3LGxMXM9Y9r4+DgURYGiKJiamlpzW4FAAC6XCy6XC/Pz8wXrEJXC4ktEDcPr9SISieRMC4fDay5b7gAKsiwDAObm5sxpTqcTDofDPMvN31YkEoGiKJAkCR6PZ811iNbDZmciqirjTuVaDCkoSRJ8Ph8CgYA5YILL5YKiKGbzsdfrNQdVyB5AIZVKYXZ21pyWSCQwOjqKYDAIp9NZMDiDoihmYc7fls/ny1kmf53sPqWJ1sLiS0QNxeFwFBQ3SZIKrrMaAydkD6Dw5JNPoq2tLWeasV6xwRmKbSs/Bl7npY1gszMREZHFWHyJiIgsxuJLRFWxurpqdwhbzurqKpqamuwOg2qAxZeIKtbW1obFxUW7w9hyFhcX0drKW3O2Iv6qRFSxrq4ufP755zlnvysrK2hurq/j+3qLqVQ8TU1N6O7utjAisgqLLxFVrKOjAx0dHTnTvvrqK+zevdumiNZWbzHVWzxknfo5BCQiIrpLVHzmK4Toy7x8StO045Vuj4jsw3wmskZFxVcI8TgySSqEOC6EUDVNSxRb/rXXXqvk44iohjaSz8xloso0VePxACGEBOCUpmnDlYdERHZiPhPVXrWu+QoAepW2RUT2Yj4T1VjJM18hxNE1Jic0Tbuct9w7AC5pmhasYnxEVEXMZ6L6UPKar6Zpk8XmCSFOAfgss4wOQK5ibERUZcxnovpQ0TVfIYQKwBjK4xleIyJqXMxnIutU5YYrIiIiKh872SAiIrKYbd1LCiG8uH1d6XFN005vdL7VMWUevzA6IHjCqg4Iyv0ehBCnrIipjN/tcWSaLq26WWcD+5K63jXPKsfUB+C4pmlPFZlv+f5dK8zlymPKW86SXC4nJuaz+ZlVz2dbznwzgSJzh6We1atOWfPtiAnAswBkYwcscteo1TEZy/XhzrU6u+MZznxHauYaoq0xZd4bd/MmMn9Mai7/7uG8mCzfv2uFuVy1mIzlLMnlDcTEfEZt8tmuZucnABg95yQA5H+BpeZbHpOmaZNZR1kqgKI/hlUxAeZNMkV7FbMynswfsWimZ6TT6/V2ZlVMADQA540jeE3TPrYgplLs2L9rhblchZgAy3O5ZEzM57Jtah+3q/hKee93bnB+LZT1mZkESVm0I5YT07pdelZZqXgOZP6lhBDvZJr3bI1J0zQdwDsAzmdiqwd27N+1wlwuT73lMsB8rpZN7eN2Fd9SzxDa8YxhuZ/ptfARjHVjEkL0rdccYnU8GZ9lEiQKoObNeSj9HXkBXNY07UDWe7ttpWdomcvlqbdcBpjP1bKpfdyu4vsv3DlaUAFc2uB8O2KCEMJrXEy36DpdqZhSQoi+zA6oWnD9o5zfzSDBmi4KS8WU3TR1EvVR9OzYv2uFuVydmKzO5XJiYj6XZ1P7uC3FN+sCfl/m/WUAEEJcWm++nTFlpp8SQkSFENFax1NOTJqmfZyZJqOw6cOOeIIApKz5Nb8TsVRMACaFEEcz85+18O5I7+3/3Dkyt3P/rhXmcnVisjqXy4yJ+ZxRi3xmJxtEREQWYycbREREFmPxJSIishiLLxERkcVYfImIiCzG4ktERGQxFl8iIiKLsfgSERFZ7P/IdV3kvrW66QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set into eval mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Initialize plots\n",
    "f, (y1_ax, y2_ax) = plt.subplots(1, 2, figsize=(8, 3))\n",
    "# Test points every 0.02 in [0,1]\n",
    "\n",
    "# Make predictions\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    test_x = torch.linspace(0, 1, 51)\n",
    "    observed_pred = likelihood(model(test_x))\n",
    "    # Get mean\n",
    "    mean = observed_pred.mean\n",
    "    # Get lower and upper confidence bounds\n",
    "    lower, upper = observed_pred.confidence_region()\n",
    "    \n",
    "# This contains predictions for both tasks, flattened out\n",
    "# The first half of the predictions is for the first task\n",
    "# The second half is for the second task\n",
    "\n",
    "# Plot training data as black stars\n",
    "y1_ax.plot(train_x.detach().numpy(), train_y[:, 0].detach().numpy(), 'k*')\n",
    "# Predictive mean as blue line\n",
    "y1_ax.plot(test_x.numpy(), mean[:, 0].numpy(), 'b')\n",
    "# Shade in confidence \n",
    "y1_ax.fill_between(test_x.numpy(), lower[:, 0].numpy(), upper[:, 0].numpy(), alpha=0.5)\n",
    "y1_ax.set_ylim([-3, 3])\n",
    "y1_ax.legend(['Observed Data', 'Mean', 'Confidence'])\n",
    "y1_ax.set_title('Observed Values (Likelihood)')\n",
    "\n",
    "# Plot training data as black stars\n",
    "y2_ax.plot(train_x.detach().numpy(), train_y[:, 1].detach().numpy(), 'k*')\n",
    "# Predictive mean as blue line\n",
    "y2_ax.plot(test_x.numpy(), mean[:, 1].numpy(), 'b')\n",
    "# Shade in confidence \n",
    "y2_ax.fill_between(test_x.numpy(), lower[:, 1].numpy(), upper[:, 1].numpy(), alpha=0.5)\n",
    "y2_ax.set_ylim([-3, 3])\n",
    "y2_ax.legend(['Observed Data', 'Mean', 'Confidence'])\n",
    "y2_ax.set_title('Observed Values (Likelihood)')\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
